module lFEMRect
	implicit none

contains

	subroutine lfemrect_0Dirichlet_derrors(ipars,rect_data,m,n,u_fun,uh,derrs,sqs)
		! get the errors in rectangles~sqs (optional) between u_fun and uh in Richardson mesh of rectangular~rect_data.
		! all the values in u_fun, uh are double precision
		! ipars, size 4: ipars(1)=0, do not give L2 norm of u_fun in derrs(1). else give L2 norm
		!				 ipars(2)=0, do not give H1-semi norm of u_fun in derrs(2); else give
		!				 ipars(3)=0, do not give L2 norm errors in derrs(3); else give
		!				 ipars(4)=0, do not give H1-semi norm errors in derrs(4); else give
		! derrs, size same as ipars: derrs(1)=0 if ipars(1)=0, else is equal to the L2 norm of u_fun
		!							 derrs(2)=0 if ipars(2)=0, else is equal to the H1-semi norm of u_fun
		!							 derrs(3)=0 if ipars(3)=0, else is equal to the L2 norm errors
		!							 derrs(4)=0 if ipars(4)=0, else give the H1-semi norm errors
		! sqs=(/ m1,m2,n1,n2 /): should 1 .le. m1 .le. m2 .le. m; 1 .le. n1 .le. n2 .le. n. means that get the errors in the union of squares~(i,j) for m1.le.i.le.m2 and n1.le.j.le.n2 
		implicit none
		interface
			subroutine u_fun(x,y,v,vx,vy)
				real(kind=8), dimension(:,:), intent(in) :: x,y
				real(kind=8), dimension(:,:), intent(out) :: v, vx, vy
			end subroutine u_fun
		end interface
		integer, intent(in) :: ipars(4), m, n
		real(kind=8), intent(in) :: rect_data(4), uh(:)
		real(kind=8), intent(out) :: derrs(4)
		integer, intent(in), optional :: sqs(4)

		real(kind=8), dimension(:,:), allocatable :: x,y, u, ux, uy
		real(kind=8), dimension(:), allocatable :: dtmp
		real(kind=8) :: ddum(1,1), O(2), width, height, hx, hy
		real(kind=8) :: gauss3hx(3), gauss3hy(3), weight3(3), basesInGs(9,4), weights9(9)
		integer :: ii, jj

		O = rect_data(1:2)	! left bottum vertice
		width = rect_data(3)
		height = rect_data(4)
		hx = width/m
		hy = height/n 
		gauss3hx = (/1.d0-sqrt(0.6),1.d0,1.d0+sqrt(0.6)/) / 2.0 ! 3 Gaussian points in [0,1]
		forall(ii=1:3,jj=1:3)
			basesInGs(ii+(jj-1)*3,1) = (1.d0- gauss3hx(ii))*(1.d0-gauss3hx(jj))
			basesInGs(ii+(jj-1)*3,2) = gauss3hx(ii)*(1.d0-gauss3hx(jj))
			basesInGs(ii+(jj-1)*3,3) = gauss3hx(ii)*gauss3hx(jj)
			basesInGs(ii+(jj-1)*3,4) = (1.d0- gauss3hx(ii))*gauss3hx(jj)
		end forall
		gauss3hx = (/1.d0-sqrt(0.6),1.d0,1.d0+sqrt(0.6)/)*(hx/2.d0) ! 3 Gaussian points in [0,hx]
		gauss3hy = (/1.d0-sqrt(0.6),1.d0,1.d0+sqrt(0.6)/)*(hy/2.d0) ! 3 Gaussian points in [0,hy]
		weight3 = (/ 5.d0/9.d0, 8.d0/9.d0, 5.d0/9.d0 /) / 2.d0	! weights
		forall(ii=1:3,jj=1:3)
			weights9(ii+(jj-1)*3) = weight3(ii)*weight3(jj)
		end forall
		allocate( dtmp(m*n) )

		derrs = 0
		if ( ipars(1).eq.0 .and. ipars(2).eq.0 .and. ipars(3).eq.0 .and. ipars(4).eq.0 ) return 
		if ( .not. present(sqs) ) then
			! get the L2 norm of u_fun
			if ( ipars(1).ne.0 ) then
				allocate(x(9,m*n))
				allocate(y(9,m*n))
				forall(ii=1:m,jj=1:n)	! value x,y	
					x(:,ii+(jj-1)*m) = (/ (O(1)+(ii-1)*hx)+gauss3hx, (O(1)+(ii-1)*hx)+gauss3hx, (O(1)+(ii-1)*hx)+gauss3hx /)
					y(1:3,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(1) 
					y(4:6,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(2)
					y(7:9,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(3)
				end forall
				allocate(u(9,m*n))
				call u_fun(x,y,u,ddum,ddum)
				call dgemv('t',9,m*n,hx*hy,u**2,9,weights9,1,0.d0,dtmp,1)
				derrs(1) = sqrt(sum(dtmp))	! the L2-norm of u_fun 
			end if
			! get L2 norm error of uh 
			if ( ipars(3).ne.0 ) then
				if ( .not.allocated(u) ) then
					allocate(x(9,m*n))
					allocate(y(9,m*n))
					forall(ii=1:m,jj=1:n)	! value x,y	
						x(:,ii+(jj-1)*m) = (/ (O(1)+(ii-1)*hx)+gauss3hx, (O(1)+(ii-1)*hx)+	gauss3hx, (O(1)+(ii-1)*hx)+gauss3hx /)
						y(1:3,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(1) 
						y(4:6,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(2)
						y(7:9,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(3)
					end forall
					allocate(u(9,m*n))
					call u_fun(x,y,u,ddum,ddum)
				end if 
				forall(ii=1:m,jj=1:n)
					u(:,ii+(jj-1)*m) = u(:,ii+(jj-1)*m) &
						& - uh(ii+(jj-1)*(m+1))*basesInGs(:,1) &
						& - uh(ii+1+(jj-1)*(m+1))*basesInGs(:,2) &
						& - uh(ii+1+jj*(m+1))*basesInGs(:,3) &
						& - uh(ii+jj*(m+1))*basesInGs(:,4)
				end forall
				call dgemv('t',9,m*n,hx*hy,u**2,9,weights9,1,0.d0,dtmp,1)
				derrs(3) = sqrt(sum(dtmp))	! the H1-semi norm of uh 
			end if
			! get the square of H1-semi norm of (u_fun)x
			if ( ipars(2).ne.0 ) then
				if ( .not.allocated(x) ) then
					allocate(x(9,m*n))
					allocate(y(9,m*n))
					forall(ii=1:m,jj=1:n)	! value x,y	
						x(:,ii+(jj-1)*m) = (/ (O(1)+(ii-1)*hx)+gauss3hx, (O(1)+(ii-1)*hx)+gauss3hx, (O(1)+(ii-1)*hx)+gauss3hx /)
						y(1:3,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(1) 
						y(4:6,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(2)
						y(7:9,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(3)
					end forall
				end if 
				if ( .not.allocated(u) )  allocate(u(9,m*n))
				call u_fun(x,y,ddum,u,ddum)
				call dgemv('t',9,m*n,hx*hy,u**2,9,weights9,1,0.d0,dtmp,1)
				derrs(2) = sum(dtmp)	! the square of H1-semi norm of (u_fun)x
			end if
			! get the square of H1-semi norm error of (u_fun)x
			if ( ipars(4).ne.0 ) then
				if ( ipars(2).eq.0 ) then
					if ( .not.allocated(x) ) then
						allocate(x(9,m*n))
						allocate(y(9,m*n))
						forall(ii=1:m,jj=1:n)	! value x,y	
							x(:,ii+(jj-1)*m) = (/ (O(1)+(ii-1)*hx)+gauss3hx, (O(1)+(ii-1)*hx)+gauss3hx, (O(1)+(ii-1)*hx)+gauss3hx /)
							y(1:3,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(1) 
							y(4:6,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(2)
							y(7:9,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(3)
						end forall
					end if 
					if ( .not.allocated(u) )  allocate(u(9,m*n))
					call u_fun(x,y,ddum,u,ddum)
				end if 
				gauss3hx = (/1.d0-sqrt(0.6),1.d0,1.d0+sqrt(0.6)/) / 2.0 ! 3 Gaussian points in [0,1]
				forall(ii=1:3,jj=1:3)
					basesInGs(ii+(jj-1)*3,1) = -(1.d0-gauss3hx(jj))/hx
					basesInGs(ii+(jj-1)*3,2) = (1.d0-gauss3hx(jj))/hx
					basesInGs(ii+(jj-1)*3,3) = gauss3hx(jj)/hx
					basesInGs(ii+(jj-1)*3,4) = -gauss3hx(jj)/hx
				end forall
				forall(ii=1:m,jj=1:n)
					u(:,ii+(jj-1)*m) = u(:,ii+(jj-1)*m) &
						& - uh(ii+(jj-1)*(m+1))*basesInGs(:,1) &
						& - uh(ii+1+(jj-1)*(m+1))*basesInGs(:,2) &
						& - uh(ii+1+jj*(m+1))*basesInGs(:,3) &
						& - uh(ii+jj*(m+1))*basesInGs(:,4)
				end forall
				call dgemv('t',9,m*n,hx*hy,u**2,9,weights9,1,0.d0,dtmp,1)
				derrs(4) = sum(dtmp)	! the square of H1-semi norm error of (u_fun)x
			end if
			! get the H1-semi norm of u_fun
			if ( ipars(2).ne.0 ) then
				if ( .not.allocated(x) ) then
					allocate(x(9,m*n))
					allocate(y(9,m*n))
					forall(ii=1:m,jj=1:n)	! value x,y	
						x(:,ii+(jj-1)*m) = (/ (O(1)+(ii-1)*hx)+gauss3hx, (O(1)+(ii-1)*hx)+gauss3hx, (O(1)+(ii-1)*hx)+gauss3hx /)
						y(1:3,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(1) 
						y(4:6,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(2)
						y(7:9,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(3)
					end forall
				end if 
				if ( .not.allocated(u) )  allocate(u(9,m*n))
				call u_fun(x,y,ddum,ddum,u)
				call dgemv('t',9,m*n,hx*hy,u**2,9,weights9,1,0.d0,dtmp,1)
				derrs(2) = sqrt( derrs(2)+sum(dtmp)	) ! the  H1-semi norm of u_fun
			end if
			! get the H1-semi norm error of u_fun
			if ( ipars(4).ne.0 ) then
				if ( ipars(2).eq.0 ) then
					if ( .not.allocated(x) ) then
						allocate(x(9,m*n))
						allocate(y(9,m*n))
						forall(ii=1:m,jj=1:n)	! value x,y	
							x(:,ii+(jj-1)*m) = (/ (O(1)+(ii-1)*hx)+gauss3hx, (O(1)+(ii-1)*hx)+gauss3hx, (O(1)+(ii-1)*hx)+gauss3hx /)
							y(1:3,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(1) 
							y(4:6,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(2)
							y(7:9,ii+(jj-1)*m) = O(2)+(jj-1)*hy + gauss3hy(3)
						end forall
					end if 
					if ( .not.allocated(u) )  allocate(u(9,m*n))
					call u_fun(x,y,ddum,ddum,u)
				end if 
				gauss3hx = (/1.d0-sqrt(0.6),1.d0,1.d0+sqrt(0.6)/) / 2.0 ! 3 Gaussian points in [0,1]
				forall(ii=1:3,jj=1:3)
					basesInGs(ii+(jj-1)*3,1) = -(1.d0- gauss3hx(ii))/hy
					basesInGs(ii+(jj-1)*3,2) = -gauss3hx(ii)/hy
					basesInGs(ii+(jj-1)*3,3) = gauss3hx(ii)/hy
					basesInGs(ii+(jj-1)*3,4) = (1.d0- gauss3hx(ii))/hy
				end forall
				forall(ii=1:m,jj=1:n)
					u(:,ii+(jj-1)*m) = u(:,ii+(jj-1)*m) &
						& - uh(ii+(jj-1)*(m+1))*basesInGs(:,1) &
						& - uh(ii+1+(jj-1)*(m+1))*basesInGs(:,2) &
						& - uh(ii+1+jj*(m+1))*basesInGs(:,3) &
						& - uh(ii+jj*(m+1))*basesInGs(:,4)
				end forall
				call dgemv('t',9,m*n,hx*hy,u**2,9,weights9,1,0.d0,dtmp,1)
				derrs(4) = sqrt( derrs(4)+sum(dtmp) )	! the H1-semi norm error of u_fun
			end if
		else
		end if 


		if ( allocated(x) ) deallocate(x)
		if ( allocated(y) ) deallocate(y)
		if ( allocated(u) ) deallocate(u)
		if ( allocated(ux) ) deallocate(ux)
		if ( allocated(uy) ) deallocate(uy)
		if ( allocated(dtmp) ) deallocate(dtmp)
	end subroutine lfemrect_0Dirichlet_derrors

	subroutine lfemrect_0Dirichlet_dsolution(rect_data,m,n,c_fun,a_fun,f_fun,u,mtype)
		!	linear continuous FEM on Cartesian Mesh
		!	get the solution of linear FEM on Rectangular~rect_data
		!	for the problem with 0 Dirichlet boundary condition:
		!		 -\div(c1*u_x,c2*u_y) + av*u = f
		!	(a,ia,ja) stored in sysmetric CSR form
		! all the functions c_fun, a_fun and f_fun are valued by double precision
		use MySparseOperator
		implicit none
		interface
			subroutine c_fun(x,y,v1,v2)
				real(kind=8), dimension(:), intent(in) :: x,y
				real(kind=8), dimension(:), intent(out) :: v1,v2
			end subroutine c_fun
			subroutine a_fun(x,y,v)
				real(kind=8), dimension(:), intent(in) :: x,y
				real(kind=8), dimension(:), intent(out) :: v
			end subroutine a_fun
			subroutine f_fun(x,y,v)
				real(kind=8), dimension(:), intent(in) :: x,y
				real(kind=8), dimension(:), intent(out) :: v
			end subroutine f_fun
		end interface
		real(kind=8) :: rect_data(4)
		integer, intent(in) :: m,n
		real(kind=8), dimension(:), allocatable, intent(out) :: u
		integer, intent(in), optional :: mtype

		integer, dimension(:), allocatable:: ia,ja
		real(kind=8), dimension(:), allocatable :: a, Fv, ut
		integer :: msglvl, ii, jj, iter

		! get the stiffness matrix without boundary condition
		call lfemrect_0Dirichlet_dstiffmtx(rect_data,m,n,c_fun,a_fun,a,ia,ja)
		! get the right vector without boundary condition
		call lfemrect_0Dirichlet_drightvec(rect_data,m,n,f_fun, Fv, 'n')
		! get the solution without boundary condition
		msglvl = 0	! do not show the message
		allocate( ut((m-1)*(n-1)) )
		if ( present(mtype) ) then
			call mso_dpardiso(mtype,msglvl,a,ia,ja,Fv,ut)
		else
			call mso_dpardiso(-2,msglvl,a,ia,ja,Fv,ut)
		end if
		allocate(u((m+1)*(n+1)))
		forall(ii=1:m-1,jj=1:n-1)
			u(ii+1+jj*(m+1)) = ut(ii+(jj-1)*(m-1))
		end forall
		forall(ii=1:m+1)
			u(ii) = 0.d0
			u(ii+n*(m+1)) = 0.d0
		end forall
		forall(jj=1:n+1)
			u(1+(jj-1)*(m+1)) = 0.d0
			u(jj*(m+1)) = 0.d0
		end forall


		if ( allocated(a) ) deallocate(a)
		if ( allocated(ia) ) deallocate(ia)
		if ( allocated(ja) ) deallocate(ja)
		if ( allocated(Fv) ) deallocate(Fv)
		if ( allocated(ut) ) deallocate(ut)
	end subroutine lfemrect_0Dirichlet_dsolution

	subroutine lfemrect_0Dirichlet_dstiffmtx(rect_data,m,n,c_fun,a_fun,a,ia,ja,b,ib,jb)
		!	linear continuous FEM on Cartesian Mesh
		!	get the stiffness matrix of linear FEM on Rectangular~rect_data
		!	for the problem: -\div(c1*u_x,c2*u_y) + av*u = f
		!	(a,ia,ja) stored in sysmetric CSR form
		use MySparseOperator
		implicit none
		interface
			subroutine c_fun(x,y,v1,v2)
				real(kind=8), dimension(:), intent(in) :: x,y
				real(kind=8), dimension(:), intent(out) :: v1,v2
			end subroutine c_fun
			subroutine a_fun(x,y,v)
				real(kind=8), dimension(:), intent(in) :: x,y
				real(kind=8), dimension(:), intent(out) :: v
			end subroutine a_fun
		end interface
		real(kind=8) :: rect_data(4)
		integer, intent(in) :: m,n
		integer, dimension(:), allocatable, intent(out) :: ia,ja
		real(kind=8), dimension(:), allocatable, intent(out) :: a
		integer, dimension(:), allocatable, optional, intent(out) :: ib,jb
		real(kind=8), dimension(:), allocatable, optional, intent(out) :: b

		real(kind=8) :: O(2), width, height, hx,hy, Axlocal(4,4),Aylocal(4,4),Alocal(4,4)
		integer :: nnz, dofs, indofs
		integer :: ii,jj
		integer, dimension(:), allocatable :: itmp, ic,jc,id,jd
		real(kind=8), dimension(:), allocatable :: val, c,d, x,y,c1m,c2m,am


		O = rect_data(1:2)
		width = rect_data(3)
		height = rect_data(4)
		hx = width/m
		hy = height/n
		Axlocal = reshape( (/	&
					&	1.d0/3.d0,	-1.d0/3.d0,	-1.d0/6.d0,	1.d0/6.d0,	&
					&	-1.d0/3.d0,	1.d0/3.d0,	1.d0/6.d0,	-1.d0/6.d0,	&
					&	-1.d0/6.d0,	1.d0/6.d0,	1.d0/3.d0,	-1.d0/3.d0,	&
					&	1.d0/6.d0,	-1.d0/6.d0,	-1.d0/3.d0,	1.d0/3.d0	/),	&
					&	(/4,4/)	)
		Aylocal = reshape( (/	&
					&	1.d0/3.d0,	1.d0/6.d0,	-1.d0/6.d0,	-1.d0/3.d0,	&
					&	1.d0/6.d0,	1.d0/3.d0,	-1.d0/3.d0,	-1.d0/6.d0,	&
					&	-1.d0/6.d0,	-1.d0/3.d0,	1.d0/3.d0,	1.d0/6.d0,	&
					&	-1.d0/3.d0,	-1.d0/6.d0,	1.d0/6.d0,	1.d0/3.d0	/),	&
					&	(/4,4/)	)
		Alocal = reshape( (/	&
					&	1.d0/9.d0,	1.d0/18.d0,	1.d0/36.d0,	1.d0/18.d0,	&
					&	1.d0/18.d0,	1.d0/9.d0,	1.d0/18.d0,	1.d0/36.d0,	&
					&	1.d0/36.d0,	1.d0/18.d0,	1.d0/9.d0,	1.d0/18.d0,	&
					&	1.d0/18.d0,	1.d0/36.d0,	1.d0/18.d0,	1.d0/9.d0	/),	&
					&	(/4,4/)	)


		if ( present(b) .and. present(ib) .and. present(jb) ) then
			dofs = (m+1)*(n+1)
			indofs = (m-1)*(n-1)
			nnz = n*(4+5*(m-1)+3) + (2*m+1)

			! value b, ib, jb			
			if ( .not. allocated(b) ) allocate(b(indofs))
			if ( .not. allocated(ib) ) allocate(ib(dofs+1))
			if ( .not. allocated(jb) ) allocate(jb(indofs))
			b = 1.d0
			forall(ii=1:indofs) jb(ii) = ii
			ib(1:m+1) = 1
			allocate(itmp(m+1))
			itmp = (/ 1, (ii,ii=1,m) /)
			forall(ii=2:n) ib((ii-1)*(m+1)+1:ii*(m+1)) = itmp+(ii-2)*(m-1)
			ib((m+1)*n+1:dofs+1) = indofs+1
			! value a, ja, ja
			allocate(a(nnz))
			allocate(ja(nnz))
			allocate(ia(dofs+1))
			
			allocate(c(nnz))
			allocate(jc(nnz))
			allocate(ic(dofs+1))

			allocate(d(m*n))
			allocate(jd(m*n))
			allocate(id(dofs+1))

			if ( allocated(itmp) ) deallocate(itmp)
			allocate(itmp(m*n))
			forall(ii=1:n) itmp((ii-1)*m+1:ii*m) = (/ (ii-1)*(m+1)+1:(ii-1)*(m+1)+m /)

			allocate(val(m*n))
			allocate(x(m*n))
			allocate(y(m*n))
			allocate(c1m(m*n))
			allocate(c2m(m*n))
			allocate(am(m*n))

			val(1:m) = (/ (O(1)+(ii+0.5)*hx,ii=0,m-1) /)
			forall(ii=1:n) 
				x((ii-1)*m+1:ii*m) = val(1:m)
				y((ii-1)*m+1:ii*m) = O(2)+(ii-0.5)*hy
			end forall

			call c_fun(x,y,c1m,c2m)
			call a_fun(x,y,am)
			deallocate(x,y)

			val = c1m*(Axlocal(1,1)*hy/hx) + c2m*(Aylocal(1,1)*hx/hy) + am*(hx*hy*Alocal(1,1))
			call mso_dcoo2csr(itmp,itmp,val,a,ia,ja,dofs,m*n)
			val = c1m*(Axlocal(1,2)*hy/hx) + c2m*(Aylocal(1,2)*hx/hy) + am*(hx*hy*Alocal(1,2))
			call mso_dcoo2csr(itmp,itmp+1,val,d,id,jd,dofs,m*n)
			call mso_dcsradd('n',dofs,dofs, a,ja,ia, 1.d0, d,jd,id, c,jc,ic)
			val = c1m*(Axlocal(1,4)*hy/hx) + c2m*(Aylocal(1,4)*hx/hy) + am*(hx*hy*Alocal(1,4))
			call mso_dcoo2csr(itmp,itmp+(m+1),val,d,id,jd,dofs,m*n)
			call mso_dcsradd('n',dofs,dofs, c,jc,ic, 1.d0, d,jd,id, a,ja,ia)
			val = c1m*(Axlocal(1,3)*hy/hx) + c2m*(Aylocal(1,3)*hx/hy) + am*(hx*hy*Alocal(1,3))
			call mso_dcoo2csr(itmp,itmp+(m+2),val,d,id,jd,dofs,m*n)
			call mso_dcsradd('n',dofs,dofs, a,ja,ia, 1.d0, d,jd,id, c,jc,ic)

			val = c1m*(Axlocal(2,2)*hy/hx) + c2m*(Aylocal(2,2)*hx/hy) + am*(hx*hy*Alocal(2,2))
			call mso_dcoo2csr(itmp+1,itmp+1,val,d,id,jd,dofs,m*n)
			call mso_dcsradd('n',dofs,dofs, c,jc,ic, 1.d0, d,jd,id, a,ja,ia)
			val = c1m*(Axlocal(2,4)*hy/hx) + c2m*(Aylocal(2,4)*hx/hy) + am*(hx*hy*Alocal(2,4))
			call mso_dcoo2csr(itmp+1,itmp+(m+1),val,d,id,jd,dofs,m*n)
			call mso_dcsradd('n',dofs,dofs, a,ja,ia, 1.d0, d,jd,id, c,jc,ic)
			val = c1m*(Axlocal(2,3)*hy/hx) + c2m*(Aylocal(2,3)*hx/hy) + am*(hx*hy*Alocal(2,3))
			call mso_dcoo2csr(itmp+1,itmp+(m+2),val,d,id,jd,dofs,m*n)
			call mso_dcsradd('n',dofs,dofs, c,jc,ic, 1.d0, d,jd,id, a,ja,ia)

			val = c1m*(Axlocal(4,4)*hy/hx) + c2m*(Aylocal(4,4)*hx/hy) + am*(hx*hy*Alocal(4,4))
			call mso_dcoo2csr(itmp+(m+1),itmp+(m+1),val,d,id,jd,dofs,m*n)
			call mso_dcsradd('n',dofs,dofs, a,ja,ia, 1.d0, d,jd,id, c,jc,ic)
			val = c1m*(Axlocal(4,3)*hy/hx) + c2m*(Aylocal(4,3)*hx/hy) + am*(hx*hy*Alocal(4,3))
			call mso_dcoo2csr(itmp+(m+1),itmp+(m+2),val,d,id,jd,dofs,m*n)
			call mso_dcsradd('n',dofs,dofs, c,jc,ic, 1.d0, d,jd,id, a,ja,ia)

			val = c1m*(Axlocal(3,3)*hy/hx) + c2m*(Aylocal(3,3)*hx/hy) + am*(hx*hy*Alocal(3,3))
			call mso_dcoo2csr(itmp+(m+2),itmp+(m+2),val,d,id,jd,dofs,m*n)
			call mso_dcsradd('n',dofs,dofs, a,ja,ia, 1.d0, d,jd,id, c,jc,ic)

			a = c
			ia = ic
			ja = jc
		else
			if ( m.lt.3 .or. n.lt.3 ) then
				write(*,*) 'm and n should be greater than or equal to 3 in lfemrect_0Dirichlet_dstiffmtx'
				stop 1
			end if 
			dofs = (m-1)*(n-1)
			nnz = 5*(m-3)*(n-2) + 4*(n-2) + 3*(n-2) + 2*(m-2) + 1
			allocate(a(nnz))
			allocate(ja(nnz))
			allocate(ia(dofs+1))
			! value ia
			forall(ii=2:m-2,jj=1:n-2)	ia(ii+(jj-1)*(m-1)) = 5
			forall(jj=1:n-2)
				ia(1+(jj-1)*(m-1)) = 4
				ia(jj*(m-1)) = 3
			end forall
			forall(ii=1:m-2) ia((n-2)*(m-1)+ii) = 2
			ia(dofs) = 1
			indofs = 0
			do ii = 1, dofs, 1
				jj = ia(ii)
				ia(ii) = indofs + 1
				indofs = indofs + jj
			end do
			ia(dofs+1) = indofs + 1
			! value ja
			forall(ii=2:m-2,jj=1:n-2)
				ja(ia(ii+(jj-1)*(m-1))) = ii+(jj-1)*(m-1)
				ja(ia(ii+(jj-1)*(m-1))+1) = ii+(jj-1)*(m-1) + 1
				ja(ia(ii+(jj-1)*(m-1))+2) = ii+(jj-1)*(m-1) + m-2
				ja(ia(ii+(jj-1)*(m-1))+3) = ii+(jj-1)*(m-1) + m-1
				ja(ia(ii+(jj-1)*(m-1))+4) = ii+(jj-1)*(m-1) + m
			end forall
			forall(jj=1:n-2)
				ja(ia(1+(jj-1)*(m-1))) = 1+(jj-1)*(m-1)
				ja(ia(1+(jj-1)*(m-1))+1) = 1+(jj-1)*(m-1) + 1
				ja(ia(1+(jj-1)*(m-1))+2) = 1+(jj-1)*(m-1) + m-1
				ja(ia(1+(jj-1)*(m-1))+3) = 1+(jj-1)*(m-1) + m
				ja(ia(jj*(m-1))) = jj*(m-1)
				ja(ia(jj*(m-1))+1) = jj*(m-1) + m-2
				ja(ia(jj*(m-1))+2) = jj*(m-1) + m-1
			end forall
			forall(ii=1:m-2)
				ja(ia((n-2)*(m-1)+ii)) = (n-2)*(m-1)+ii
				ja(ia((n-2)*(m-1)+ii)+1) = (n-2)*(m-1)+ii + 1
			end forall
			ja(ia(dofs)) = dofs
			! value a
			allocate(x(m*n))
			allocate(y(m*n))
			allocate(c1m(m*n))
			allocate(c2m(m*n))
			allocate(am(m*n))
			forall(ii=1:m,jj=1:n) 
				x(ii+(jj-1)*m) = O(1)+(ii-0.5)*hx
				y(ii+(jj-1)*m) = O(2)+(jj-0.5)*hy
			end forall
			call c_fun(x,y,c1m,c2m)
			call a_fun(x,y,am)
			deallocate(x,y)
			forall(ii=2:m-2,jj=1:n-2)
				a(ia(ii+(jj-1)*(m-1))) = c1m(ii+(jj-1)*m)*Axlocal(3,3)*hy/hx + c2m(ii+(jj-1)*m)*Aylocal(3,3)*hx/hy + am(ii+(jj-1)*m)*hx*hy*Alocal(3,3) +&
					& c1m(ii+1+(jj-1)*m)*Axlocal(4,4)*hy/hx + c2m(ii+1+(jj-1)*m)*Aylocal(4,4)*hx/hy + am(ii+1+(jj-1)*m)*hx*hy*Alocal(4,4) +&
					& c1m(ii+1+jj*m)*Axlocal(1,1)*hy/hx + c2m(ii+1+jj*m)*Aylocal(1,1)*hx/hy + am(ii+1+jj*m)*hx*hy*Alocal(1,1) +&
					& c1m(ii+jj*m)*Axlocal(2,2)*hy/hx + c2m(ii+jj*m)*Aylocal(2,2)*hx/hy + am(ii+jj*m)*hx*hy*Alocal(2,2)
				a(ia(ii+(jj-1)*(m-1))+1) = c1m(ii+1+(jj-1)*m)*Axlocal(3,4)*hy/hx + c2m(ii+1+(jj-1)*m)*Aylocal(3,4)*hx/hy + am(ii+1+(jj-1)*m)*hx*hy*Alocal(3,4) +&
					& c1m(ii+1+jj*m)*Axlocal(1,2)*hy/hx + c2m(ii+1+jj*m)*Aylocal(1,2)*hx/hy + am(ii+1+jj*m)*hx*hy*Alocal(1,2)
				a(ia(ii+(jj-1)*(m-1))+2) = c1m(ii+jj*m)*Axlocal(2,4)*hy/hx + c2m(ii+jj*m)*Aylocal(2,4)*hx/hy + am(ii+jj*m)*hx*hy*Alocal(2,4)
				a(ia(ii+(jj-1)*(m-1))+3) = c1m(ii+jj*m)*Axlocal(2,3)*hy/hx + c2m(ii+jj*m)*Aylocal(2,3)*hx/hy + am(ii+jj*m)*hx*hy*Alocal(2,3) +&
					& c1m(ii+1+jj*m)*Axlocal(1,4)*hy/hx + c2m(ii+1+jj*m)*Aylocal(1,4)*hx/hy + am(ii+1+jj*m)*hx*hy*Alocal(1,4)
				a(ia(ii+(jj-1)*(m-1))+4) = c1m(ii+1+jj*m)*Axlocal(3,1)*hy/hx + c2m(ii+1+jj*m)*Aylocal(3,1)*hx/hy + am(ii+1+jj*m)*hx*hy*Alocal(3,1)	
			end forall
			forall(jj=1:n-2)
				a(ia(1+(jj-1)*(m-1))) = c1m(1+(jj-1)*m)*Axlocal(3,3)*hy/hx + c2m(1+(jj-1)*m)*Aylocal(3,3)*hx/hy + am(1+(jj-1)*m)*hx*hy*Alocal(3,3) +&
					& c1m(1+1+(jj-1)*m)*Axlocal(4,4)*hy/hx + c2m(1+1+(jj-1)*m)*Aylocal(4,4)*hx/hy + am(1+1+(jj-1)*m)*hx*hy*Alocal(4,4) +&
					& c1m(1+1+jj*m)*Axlocal(1,1)*hy/hx + c2m(1+1+jj*m)*Aylocal(1,1)*hx/hy + am(1+1+jj*m)*hx*hy*Alocal(1,1) +&
					& c1m(1+jj*m)*Axlocal(2,2)*hy/hx + c2m(1+jj*m)*Aylocal(2,2)*hx/hy + am(1+jj*m)*hx*hy*Alocal(2,2)
				a(ia(1+(jj-1)*(m-1))+1) = c1m(1+1+(jj-1)*m)*Axlocal(3,4)*hy/hx + c2m(1+1+(jj-1)*m)*Aylocal(3,4)*hx/hy + am(1+1+(jj-1)*m)*hx*hy*Alocal(3,4) +&
					& c1m(1+1+jj*m)*Axlocal(1,2)*hy/hx + c2m(1+1+jj*m)*Aylocal(1,2)*hx/hy + am(1+1+jj*m)*hx*hy*Alocal(1,2)
				a(ia(1+(jj-1)*(m-1))+2) = c1m(1+jj*m)*Axlocal(2,3)*hy/hx + c2m(1+jj*m)*Aylocal(2,3)*hx/hy + am(1+jj*m)*hx*hy*Alocal(2,3) +&
					& c1m(1+1+jj*m)*Axlocal(1,4)*hy/hx + c2m(1+1+jj*m)*Aylocal(1,4)*hx/hy + am(1+1+jj*m)*hx*hy*Alocal(1,4)
				a(ia(1+(jj-1)*(m-1))+3) = c1m(1+1+jj*m)*Axlocal(3,1)*hy/hx + c2m(1+1+jj*m)*Aylocal(3,1)*hx/hy + am(1+1+jj*m)*hx*hy*Alocal(3,1)

				a(ia(jj*(m-1))) = c1m(m-1+(jj-1)*m)*Axlocal(3,3)*hy/hx + c2m(m-1+(jj-1)*m)*Aylocal(3,3)*hx/hy + am(m-1+(jj-1)*m)*hx*hy*Alocal(3,3) +&
					& c1m(m-1+1+(jj-1)*m)*Axlocal(4,4)*hy/hx + c2m(m-1+1+(jj-1)*m)*Aylocal(4,4)*hx/hy + am(m-1+1+(jj-1)*m)*hx*hy*Alocal(4,4) +&
					& c1m(m-1+1+jj*m)*Axlocal(1,1)*hy/hx + c2m(m-1+1+jj*m)*Aylocal(1,1)*hx/hy + am(m-1+1+jj*m)*hx*hy*Alocal(1,1) +&
					& c1m(m-1+jj*m)*Axlocal(2,2)*hy/hx + c2m(m-1+jj*m)*Aylocal(2,2)*hx/hy + am(m-1+jj*m)*hx*hy*Alocal(2,2)
				a(ia(jj*(m-1))+1) = c1m(m-1+jj*m)*Axlocal(2,4)*hy/hx + c2m(m-1+jj*m)*Aylocal(2,4)*hx/hy + am(m-1+jj*m)*hx*hy*Alocal(2,4)
				a(ia(jj*(m-1))+2) = c1m(m-1+jj*m)*Axlocal(2,3)*hy/hx + c2m(m-1+jj*m)*Aylocal(2,3)*hx/hy + am(m-1+jj*m)*hx*hy*Alocal(2,3) +&
					& c1m(m-1+1+jj*m)*Axlocal(1,4)*hy/hx + c2m(m-1+1+jj*m)*Aylocal(1,4)*hx/hy + am(m-1+1+jj*m)*hx*hy*Alocal(1,4)
			end forall
			jj = n-1
			forall(ii=1:m-2)
				a(ia((n-2)*(m-1)+ii)) = c1m(ii+(jj-1)*m)*Axlocal(3,3)*hy/hx + c2m(ii+(jj-1)*m)*Aylocal(3,3)*hx/hy + am(ii+(jj-1)*m)*hx*hy*Alocal(3,3) +&
					& c1m(ii+1+(jj-1)*m)*Axlocal(4,4)*hy/hx + c2m(ii+1+(jj-1)*m)*Aylocal(4,4)*hx/hy + am(ii+1+(jj-1)*m)*hx*hy*Alocal(4,4) +&
					& c1m(ii+1+jj*m)*Axlocal(1,1)*hy/hx + c2m(ii+1+jj*m)*Aylocal(1,1)*hx/hy + am(ii+1+jj*m)*hx*hy*Alocal(1,1) +&
					& c1m(ii+jj*m)*Axlocal(2,2)*hy/hx + c2m(ii+jj*m)*Aylocal(2,2)*hx/hy + am(ii+jj*m)*hx*hy*Alocal(2,2)
				a(ia((n-2)*(m-1)+ii)+1) = c1m(ii+1+(jj-1)*m)*Axlocal(3,4)*hy/hx + c2m(ii+1+(jj-1)*m)*Aylocal(3,4)*hx/hy + am(ii+1+(jj-1)*m)*hx*hy*Alocal(3,4) +&
					& c1m(ii+1+jj*m)*Axlocal(1,2)*hy/hx + c2m(ii+1+jj*m)*Aylocal(1,2)*hx/hy + am(ii+1+jj*m)*hx*hy*Alocal(1,2)
			end forall
			ii = m-1
			jj = n-1
			a(ia(dofs)) = c1m(ii+(jj-1)*m)*Axlocal(3,3)*hy/hx + c2m(ii+(jj-1)*m)*Aylocal(3,3)*hx/hy + am(ii+(jj-1)*m)*hx*hy*Alocal(3,3) +&
					& c1m(ii+1+(jj-1)*m)*Axlocal(4,4)*hy/hx + c2m(ii+1+(jj-1)*m)*Aylocal(4,4)*hx/hy + am(ii+1+(jj-1)*m)*hx*hy*Alocal(4,4) +&
					& c1m(ii+1+jj*m)*Axlocal(1,1)*hy/hx + c2m(ii+1+jj*m)*Aylocal(1,1)*hx/hy + am(ii+1+jj*m)*hx*hy*Alocal(1,1) +&
					& c1m(ii+jj*m)*Axlocal(2,2)*hy/hx + c2m(ii+jj*m)*Aylocal(2,2)*hx/hy + am(ii+jj*m)*hx*hy*Alocal(2,2)
		end if

		if ( allocated(x) ) deallocate(x)
		if ( allocated(y) ) deallocate(y)
		if ( allocated(c1m) ) deallocate(c1m)
		if ( allocated(c2m) ) deallocate(c2m)
		if ( allocated(am) ) deallocate(am)
		if ( allocated(itmp) ) deallocate(itmp)
		if ( allocated(val) ) deallocate(val)
		if ( allocated(c) ) deallocate(c)
		if ( allocated(ic) ) deallocate(ic)
		if ( allocated(jc) ) deallocate(jc)
		if ( allocated(d) ) deallocate(d)
		if ( allocated(id) ) deallocate(id)
		if ( allocated(jd) ) deallocate(jd)      
	end subroutine lfemrect_0Dirichlet_dstiffmtx


	subroutine lfemrect_0Dirichlet_drightvec(rect_data,m,n,f_fun, Fv, str)
		! linear continuous FEM on Cartesian Mesh
		! get the right vector of linear FEM on Rectangular~rect_data
		! for the problem: -\div(c1*u_x,c2*u_y) + a*u = f
		implicit none
		interface
			subroutine f_fun(x,y,v)
				real(kind=8), dimension(:), intent(in) :: x,y
				real(kind=8), dimension(:), intent(out) :: v
			end subroutine f_fun
		end interface
		real(kind=8) :: rect_data(4)
		integer, intent(in) :: m,n
		real(kind=8), dimension(:), allocatable :: Fv
		character, optional :: str(4)

		real(kind=8) :: O(2), width, height, hx, hy
		real(kind=8), dimension(:), allocatable :: x,y, fm
		integer, dimension(:), allocatable :: itmp
		integer :: ii,jj

		O = rect_data(1:2)
		width = rect_data(3)
		height = rect_data(4)
		hx = width/m
		hy = height/n

		if ( (.not.present(str)) ) then
			allocate(x(m*n))
			allocate(y(m*n))
			allocate(fm(m*n))
			fm(1:m) = (/ (O(1)+(ii+0.5)*hx,ii=0,m-1) /)
			forall(ii=1:n) 
				x((ii-1)*m+1:ii*m) = fm(1:m)
				y((ii-1)*m+1:ii*m) = O(2)+(ii-0.5)*hy
			end forall
			call f_fun(x,y,fm)
			deallocate(x,y)

			if ( allocated(itmp) ) deallocate(itmp)
			allocate(itmp(m*n))
			forall(ii=1:n) itmp((ii-1)*m+1:ii*m) = (/ (ii-1)*(m+1)+1:(ii-1)*(m+1)+m /)

			allocate(Fv((m+1)*(n+1)))
			Fv = 0
			Fv(itmp) = Fv(itmp) + fm*(hx*hy/4)
			Fv(itmp+1) = Fv(itmp+1) + fm*(hx*hy/4)
			Fv(itmp+(m+1)) = Fv(itmp+(m+1)) + fm*(hx*hy/4)
			Fv(itmp+(m+2)) = Fv(itmp+(m+2)) + fm*(hx*hy/4)
		else
			allocate(x(m*n))
			allocate(y(m*n))
			allocate(fm(m*n))
			fm(1:m) = (/ (O(1)+(ii+0.5)*hx,ii=0,m-1) /)
			forall(ii=1:n) 
				x((ii-1)*m+1:ii*m) = fm(1:m)
				y((ii-1)*m+1:ii*m) = O(2)+(ii-0.5)*hy
			end forall
			call f_fun(x,y,fm)
			deallocate(x,y)
			allocate(Fv((m-1)*(n-1)))
			forall(ii=1:m-1,jj=1:n-1)
				Fv(ii+(jj-1)*(m-1)) = ( fm(ii+(jj-1)*m) + fm(ii+(jj-1)*m+1) + fm(ii+jj*m) + fm(ii+1+jj*m) ) * hx*hy/4.0
			end forall
		end if

		if ( allocated(x) ) deallocate(x)
		if ( allocated(y) ) deallocate(y)
		if ( allocated(fm) ) deallocate(fm)  
		if ( allocated(itmp) ) deallocate(itmp) 
	end subroutine lfemrect_0Dirichlet_drightvec




end module lFEMRect

module FEMTri

contains

	subroutine femtrg(p,t,ar,gx,gy,sl,nx,ny)
		use functions_matlab
		implicit none
		real(kind=8), dimension(:,:), intent(in) :: p
		integer, dimension(:,:), intent(in) :: t
		real(kind=8), dimension(:), intent(out) :: ar
		real(kind=8), dimension(:,:), intent(out), optional :: gx, gy, nx, ny, sl

		integer :: st(2)
		real(kind=8), dimension(:,:), allocatable :: dx, dy

		st = shape(t)
		allocate(dx(3,st(2)))
		allocate(dy(3,st(2)))

		dx(1,:) = p(1,t(2,:)) - p(1,t(3,:))
		dx(2,:) = p(1,t(3,:)) - p(1,t(1,:))
		dx(3,:) = p(1,t(1,:)) - p(1,t(2,:))
		dy(1,:) = p(2,t(2,:)) - p(2,t(3,:))
		dy(2,:) = p(2,t(3,:)) - p(2,t(1,:))
		dy(3,:) = p(2,t(1,:)) - p(2,t(2,:))

		ar = 0.5*abs(-dx(1,:)*dy(2,:)+dx(2,:)*dy(1,:))
		
		if ( present(gx) ) then
			call bsxfun_rmv_divide( dy, 2*ar, gx, 'R' )
			if ( present(gy) ) then
				call bsxfun_rmv_divide( dx, 2*ar, gy, 'R' )
				if ( present(sl) ) then
					sl = sqrt(dx**2+dy**2)
					if ( present(nx) ) then
						nx = - dy/sl
						if ( present(ny) ) then
							ny = dx/sl
						end if 
					end if 
				end if 
			end if 
		end if 
		deallocate(dx,dy)
	end subroutine femtrg

	subroutine RegTri_unit_square(n,p,e,t)
		! ragular triangulation of unit square [0,1]x[0,1]
		integer, intent(in) :: n
		integer, dimension(:,:), intent(out) :: e, t
		real(kind=8), dimension(:,:), intent(out) :: p 

		integer :: i
		real(kind=8) :: h
		real(kind=8), dimension(:), allocatable :: tmp
		integer, dimension(:), allocatable :: itmp

		h = 1.0/n
		allocate( tmp(n+1) )
		tmp = (/ ((i-1)*h, i=1,n+1) /)
		forall (i=1:n+1)
			p( 1, (i-1)*(n+1)+1:i*(n+1) ) = tmp
			p( 2, (i-1)*(n+1)+1:i*(n+1) ) = (i-1)*h	
		end forall
		deallocate(tmp)
	
		allocate(itmp(n**2))
		forall (i=1:n)
			itmp(1+(i-1)*n:i*n) = (/ (i-1)*(n+1)+1:(i-1)*(n+1)+n /) 
		end forall
		i = 2*n**2
		t(1,1:i:2) = itmp
		t(2,1:i:2) = itmp + (n+2)
		t(3,1:i:2) = itmp + (n+1)
		t(1,2:i:2) = itmp
		t(2,2:i:2) = itmp + 1
		t(3,2:i:2) = itmp + (n+2)
		deallocate(itmp)
	end subroutine RegTri_unit_square
end module FEMTri